module EDSurfaceRadiationMod
   
   !-------------------------------------------------------------------------------------
   ! EDSurfaceRadiation
   !
   ! This module contains function and type definitions for all things related
   ! to radiative transfer in ED modules at the land surface.
   !
   !-------------------------------------------------------------------------------------

#include "shr_assert.h"
   
  use EDTypesMod        , only : ed_patch_type, ed_site_type
  use EDTypesMod        , only : maxPatchesPerSite
  use EDTypesMod        , only : maxpft
  use FatesConstantsMod , only : r8 => fates_r8
  use FatesConstantsMod , only : itrue
  use FatesConstantsMod , only : pi_const
  use FatesInterfaceMod , only : bc_in_type
  use FatesInterfaceMod , only : bc_out_type
  use FatesInterfaceMod , only : hlm_numSWb
  use FatesInterfaceMod , only : numpft
  use EDTypesMod        , only : maxSWb
  use EDTypesMod        , only : nclmax
  use EDTypesMod        , only : nlevleaf
  use EDTypesMod        , only : n_rad_stream_types
  use EDTypesMod        , only : idiffuse
  use EDTypesMod        , only : idirect
  use EDTypesMod        , only : ivis
  use EDTypesMod        , only : inir
  use EDTypesMod        , only : ipar
  use EDCanopyStructureMod, only: calc_areaindex
  use FatesGlobals      , only : fates_log
  use FatesGlobals, only      : endrun => fates_endrun

  ! CIME globals
  use shr_log_mod       , only : errMsg => shr_log_errMsg

  implicit none

  private
  public :: ED_Norman_Radiation  ! Surface albedo and two-stream fluxes
  public :: PatchNormanRadiation
  public :: ED_SunShadeFracs
  
  logical :: debug = .false.  ! for debugging this module

  
  real(r8), public  :: albice(maxSWb) = &       ! albedo land ice by waveband (1=vis, 2=nir)
        (/ 0.80_r8, 0.55_r8 /)

contains
   
  subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
      !

      !
      ! !USES:
      use EDPftvarcon       , only : EDPftvarcon_inst
      use EDtypesMod        , only : ed_patch_type
      use EDTypesMod        , only : ed_site_type


      ! !ARGUMENTS:
      
      integer,            intent(in)            :: nsites 
      type(ed_site_type), intent(inout), target :: sites(nsites)      ! FATES site vector
      type(bc_in_type),   intent(in)            :: bc_in(nsites)
      type(bc_out_type),  intent(inout)         :: bc_out(nsites)


      ! !LOCAL VARIABLES:
      integer :: s                                   ! site loop counter
      integer :: ifp                                 ! patch loop counter
      integer :: ib                                  ! radiation broad band counter
      type(ed_patch_type), pointer :: currentPatch   ! patch pointer

      !-----------------------------------------------------------------------
      ! -------------------------------------------------------------------------------
      ! TODO (mv, 2014-10-29) the filter here is different than below 
      ! this is needed to have the VOC's be bfb - this needs to be
      ! re-examined int he future
      ! RGK,2016-08-06: FATES is still incompatible with VOC emission module
      ! -------------------------------------------------------------------------------


        do s = 1, nsites

           ifp = 0
           currentpatch => sites(s)%oldest_patch
           do while (associated(currentpatch))  
              ifp = ifp+1
              
              currentPatch%f_sun      (:,:,:) = 0._r8
              currentPatch%fabd_sun_z (:,:,:) = 0._r8
              currentPatch%fabd_sha_z (:,:,:) = 0._r8
              currentPatch%fabi_sun_z (:,:,:) = 0._r8
              currentPatch%fabi_sha_z (:,:,:) = 0._r8
              currentPatch%fabd       (:)     = 0._r8
              currentPatch%fabi       (:)     = 0._r8

              ! zero diagnostic radiation profiles
              currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8
              currentPatch%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0._r8
              currentPatch%nrmlzd_parprof_dir_z(:,:,:) = 0._r8
              currentPatch%nrmlzd_parprof_dif_z(:,:,:) = 0._r8

              currentPatch%solar_zenith_flag         = bc_in(s)%filter_vegzen_pa(ifp)
              currentPatch%solar_zenith_angle        = bc_in(s)%coszen_pa(ifp)
              currentPatch%gnd_alb_dif(1:hlm_numSWb) = bc_in(s)%albgr_dif_rb(1:hlm_numSWb)
              currentPatch%gnd_alb_dir(1:hlm_numSWb) = bc_in(s)%albgr_dir_rb(1:hlm_numSWb)

              if(currentPatch%solar_zenith_flag )then

                 bc_out(s)%albd_parb(ifp,:)            = 0._r8  ! output HLM
                 bc_out(s)%albi_parb(ifp,:)            = 0._r8  ! output HLM
                 bc_out(s)%fabi_parb(ifp,:)            = 0._r8  ! output HLM
                 bc_out(s)%fabd_parb(ifp,:)            = 0._r8  ! output HLM
                 bc_out(s)%ftdd_parb(ifp,:)            = 1._r8 ! output HLM
                 bc_out(s)%ftid_parb(ifp,:)            = 1._r8 ! output HLM
                 bc_out(s)%ftii_parb(ifp,:)            = 1._r8 ! output HLM
                 
                 if (maxval(currentPatch%nrad(1,:))==0)then
                    !there are no leaf layers in this patch. it is effectively bare ground. 
                    ! no radiation is absorbed  
                    bc_out(s)%fabd_parb(ifp,:) = 0.0_r8
                    bc_out(s)%fabi_parb(ifp,:) = 0.0_r8
                    do ib = 1,hlm_numSWb
                       bc_out(s)%albd_parb(ifp,ib) = bc_in(s)%albgr_dir_rb(ib)
                       bc_out(s)%albi_parb(ifp,ib) = bc_in(s)%albgr_dif_rb(ib)
                       bc_out(s)%ftdd_parb(ifp,ib)= 1.0_r8
                       bc_out(s)%ftid_parb(ifp,ib)= 1.0_r8
                       bc_out(s)%ftii_parb(ifp,ib)= 1.0_r8
                    enddo

                 else

                    call PatchNormanRadiation (currentPatch, &
                         bc_out(s)%albd_parb(ifp,:), &
                         bc_out(s)%albi_parb(ifp,:), &
                         bc_out(s)%fabd_parb(ifp,:), &
                         bc_out(s)%fabi_parb(ifp,:), &
                         bc_out(s)%ftdd_parb(ifp,:), &
                         bc_out(s)%ftid_parb(ifp,:), &
                         bc_out(s)%ftii_parb(ifp,:))
                    
                    
                 endif ! is there vegetation? 
                 
              end if    ! if the vegetation and zenith filter is active
              currentPatch => currentPatch%younger
           end do       ! Loop linked-list patches
        enddo           ! Loop Sites
        
      return
    end subroutine ED_Norman_Radiation
    

 ! ======================================================================================

  subroutine PatchNormanRadiation (currentPatch, &
                                   albd_parb_out, &   ! (ifp,ib)
                                   albi_parb_out, &   ! (ifp,ib)
                                   fabd_parb_out, &   ! (ifp,ib)
                                   fabi_parb_out, &   ! (ifp,ib)
                                   ftdd_parb_out, &   ! (ifp,ib)
                                   ftid_parb_out, &   ! (ifp,ib)
                                   ftii_parb_out)     ! (ifp,ib)

    ! -----------------------------------------------------------------------------------
    !
    ! This routine performs the Norman Radiation scattering for each patch.
    !
    ! -----------------------------------------------------------------------------------

    !
    ! !USES:
    use EDPftvarcon       , only : EDPftvarcon_inst
    use EDtypesMod        , only : ed_patch_type

    ! -----------------------------------------------------------------------------------
    ! !ARGUMENTS:
    ! -----------------------------------------------------------------------------------
    
    type(ed_patch_type), intent(inout), target :: currentPatch
    real(r8), intent(inout) :: albd_parb_out(hlm_numSWb)
    real(r8), intent(inout) :: albi_parb_out(hlm_numSWb)
    real(r8), intent(inout) :: fabd_parb_out(hlm_numSWb)
    real(r8), intent(inout) :: fabi_parb_out(hlm_numSWb)
    real(r8), intent(inout) :: ftdd_parb_out(hlm_numSWb)
    real(r8), intent(inout) :: ftid_parb_out(hlm_numSWb)
    real(r8), intent(inout) :: ftii_parb_out(hlm_numSWb)

    ! Locals
    ! -----------------------------------------------------------------------------------

    integer  :: radtype, L, ft, j
    integer  :: iter                                          ! Iteration index
    integer  :: irep                                          ! Flag to exit iteration loop
    real(r8) :: sb
    real(r8) :: error                                         ! Error check
    real(r8) :: down_rad, up_rad                              ! Iterative solution do Dif_dn and Dif_up
    real(r8) :: ftweight(nclmax,maxpft,nlevleaf)
    real(r8) :: k_dir(maxpft)                              ! Direct beam extinction coefficient
    real(r8) :: tr_dir_z(nclmax,maxpft,nlevleaf)         ! Exponential transmittance of direct beam radiation through a single layer
    real(r8) :: tr_dif_z(nclmax,maxpft,nlevleaf)         ! Exponential transmittance of diffuse radiation through a single layer
    real(r8) :: weighted_dir_tr(nclmax)
    real(r8) :: weighted_fsun(nclmax)
    real(r8) :: weighted_dif_ratio(nclmax,maxSWb)
    real(r8) :: weighted_dif_down(nclmax)
    real(r8) :: weighted_dif_up(nclmax)
    real(r8) :: refl_dif(nclmax,maxpft,nlevleaf,maxSWb)  ! Term for diffuse radiation reflected by laye
    real(r8) :: tran_dif(nclmax,maxpft,nlevleaf,maxSWb)  ! Term for diffuse radiation transmitted by layer
    real(r8) :: dif_ratio(nclmax,maxpft,nlevleaf,maxSWb) ! Ratio of upward to forward diffuse fluxes
    real(r8) :: Dif_dn(nclmax,maxpft,nlevleaf)           ! Forward diffuse flux onto canopy layer J (W/m**2 ground area)
    real(r8) :: Dif_up(nclmax,maxpft,nlevleaf)           ! Upward diffuse flux above canopy layer J (W/m**2 ground area)
    real(r8) :: lai_change(nclmax,maxpft,nlevleaf)       ! Forward diffuse flux onto canopy layer J (W/m**2 ground area)
    real(r8) :: f_not_abs(maxpft,maxSWb)                   ! Fraction reflected + transmitted. 1-absorbtion.
    real(r8) :: Abs_dir_z(maxpft,nlevleaf)
    real(r8) :: Abs_dif_z(maxpft,nlevleaf)
    real(r8) :: abs_rad(maxSWb)                               !radiation absorbed by soil
    real(r8) :: tr_soili                                      ! Radiation transmitted to the soil surface.
    real(r8) :: tr_soild                                      ! Radiation transmitted to the soil surface.
    real(r8) :: phi1b(maxpft)                                 ! Radiation transmitted to the soil surface.
    real(r8) :: phi2b(maxpft)
    real(r8) :: laisum                                        ! cumulative lai+sai for canopy layer (at middle of layer)
    real(r8) :: angle
    
    real(r8),parameter :: tolerance = 0.000000001_r8
    
    
    integer, parameter :: max_diag_nlevleaf = 4
    integer, parameter :: diag_nlevleaf = min(nlevleaf,max_diag_nlevleaf)  ! for diagnostics, write a small number of leaf layers
    
    real(r8) :: denom
    real(r8) :: lai_reduction(nclmax)
    
    integer  :: fp,iv,s      ! array indices
    integer  :: ib               ! waveband number
    real(r8) :: cosz             ! 0.001 <= coszen <= 1.000
    real(r8) :: chil
    real(r8) :: gdir
    

    real(r8), parameter :: forc_dir(n_rad_stream_types) = (/ 1.0_r8, 0.0_r8 /)   ! These are binary switches used
    real(r8), parameter :: forc_dif(n_rad_stream_types) = (/ 0.0_r8, 1.0_r8 /)   ! to turn off and on radiation streams


    
    associate(&
         rhol         =>    EDPftvarcon_inst%rhol                     , & ! Input:  [real(r8) (:)   ] leaf reflectance: 1=vis, 2=nir
         rhos         =>    EDPftvarcon_inst%rhos                     , & ! Input:  [real(r8) (:)   ] stem reflectance: 1=vis, 2=nir
         taul         =>    EDPftvarcon_inst%taul                     , & ! Input:  [real(r8) (:)   ] leaf transmittance: 1=vis, 2=nir
         taus         =>    EDPftvarcon_inst%taus                     , & ! Input:  [real(r8) (:)   ] stem transmittance: 1=vis, 2=nir
         xl           =>    EDPftvarcon_inst%xl                       , & ! Input:  [real(r8) (:)   ] ecophys const - leaf/stem orientation index
         clumping_index  => EDPftvarcon_inst%clumping_index) 


      ! Initialize local arrays

      weighted_dir_tr(:)   = 0._r8
      weighted_dif_down(:) = 0._r8
      weighted_dif_up(:)   = 0._r8

      tr_dir_z(:,:,:)      = 0._r8   
      tr_dif_z(:,:,:)      = 0._r8
      lai_change(:,:,:)    = 0._r8
      Dif_up(:,:,:)        = 0._r8
      Dif_dn(:,:,:)        = 0._r8
      refl_dif(:,:,:,:)    = 0.0_r8
      tran_dif(:,:,:,:)    = 0.0_r8
      dif_ratio(:,:,:,:)   = 0.0_r8


      ! Initialize the ouput arrays
      ! ---------------------------------------------------------------------------------
      albd_parb_out(1:hlm_numSWb) = 0.0_r8
      albi_parb_out(1:hlm_numSWb) = 0.0_r8
      fabd_parb_out(1:hlm_numSWb) = 0.0_r8
      fabi_parb_out(1:hlm_numSWb) = 0.0_r8
      ftdd_parb_out(1:hlm_numSWb) = 1.0_r8
      ftid_parb_out(1:hlm_numSWb) = 1.0_r8
      ftii_parb_out(1:hlm_numSWb) = 1.0_r8
                 
      ! Is this pft/canopy layer combination present in this patch?
      
      do L = 1,nclmax
         do ft = 1,numpft
            currentPatch%canopy_mask(L,ft) = 0
            do  iv = 1, currentPatch%nrad(L,ft)
               if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then
                  currentPatch%canopy_mask(L,ft) = 1
                  !I think 'present' is only used here...
               endif
            end do !iv
         end do !ft
      end do !L


      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
      ! Direct beam extinction coefficient, k_dir. PFT specific.
      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
      cosz = max(0.001_r8, currentPatch%solar_zenith_angle ) !copied from previous radiation code...
      do ft = 1,numpft
         sb = (90._r8 - (acos(cosz)*180._r8/pi_const)) * (pi_const / 180._r8)
         chil = xl(ft) !min(max(xl(ft), -0.4_r8), 0.6_r8 )
         if ( abs(chil) <= 0.01_r8) then
            chil  = 0.01_r8
         end if
         phi1b(ft) = 0.5_r8 - 0.633_r8*chil - 0.330_r8*chil*chil
         phi2b(ft) = 0.877_r8 * (1._r8 - 2._r8*phi1b(ft)) !0 = horiz leaves, 1 - vert leaves.
         gdir = phi1b(ft) + phi2b(ft) * sin(sb)
         !how much direct light penetrates a singleunit of lai?
         k_dir(ft) = clumping_index(ft) * gdir / sin(sb)
      end do !FT
         



      !do this once for one unit of diffuse, and once for one unit of direct radiation
      do radtype = 1, n_rad_stream_types 

         ! Extract information that needs to be provided by ED into local array.
         ! RGK: NOT SURE WHY WE NEED FTWEIGHT ...
         ! ------------------------------------------------------------------------------

         ftweight(:,:,:) = 0._r8
         do L = 1,currentPatch%NCL_p
            do ft = 1,numpft
               do  iv = 1, currentPatch%nrad(L,ft)
                  !this is already corrected for area in CLAP
                  ftweight(L,ft,iv) = currentPatch%canopy_area_profile(L,ft,iv) 
               end do  !iv
            end do  !ft1
         end do  !L
         if (sum(ftweight(1,:,1))<0.999_r8)then
            write(fates_log(),*) 'canopy not full',ftweight(1,:,1)
         endif
         if (sum(ftweight(1,:,1))>1.0001_r8)then
            write(fates_log(),*) 'canopy too full',ftweight(1,:,1)
         endif
              
         do L = 1,currentPatch%NCL_p !start at the top canopy layer (1 is the top layer.)

            weighted_dir_tr(L)                 = 0.0_r8
            weighted_fsun(L)                   = 0._r8
            weighted_dif_ratio(L,1:hlm_numSWb) = 0._r8
            
            !Each canopy layer (canopy, understorey) has multiple 'parallel' pft's

            do ft =1,numpft

               if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves.
                  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                  ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z.
                  ! Estimated do nine sky angles in increments of 10 degrees
                  ! PFT specific...
                  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                  tr_dif_z(L,ft,:) = 0._r8
                  do iv = 1,currentPatch%nrad(L,ft)
                     do j = 1,9
                        angle = (5._r8 + real(j - 1,r8) * 10._r8) * pi_const / 180._r8
                        gdir = phi1b(ft) + phi2b(ft) * sin(angle)
                        tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) + exp(-clumping_index(ft) * &
                             gdir / sin(angle) * &
                             (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv))) * &
                             sin(angle)*cos(angle)
                     end do
                     
                     tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) * 2._r8 * (10._r8 * pi_const / 180._r8)
                     
                  end do
                  

                  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                  ! Direct beam transmittance, tr_dir_z, uses cumulative LAI above layer J to give
                  ! unscattered direct beam onto layer J. do each PFT section.
                  ! This is just an  decay curve based on k_dir. (leaf & sun angle)
                  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                  if (L==1)then
                     tr_dir_z(L,ft,1) = 1._r8
                  else
                     tr_dir_z(L,ft,1)  = weighted_dir_tr(L-1)
                  endif
                  laisum = 0.00_r8                         
                  !total direct beam getting to the bottom of the top canopy.
                  do iv = 1,currentPatch%nrad(L,ft)
                     laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv)
                     lai_change(L,ft,iv) = 0.0_r8
                     if (( ftweight(L,ft,iv+1)  >  0.0_r8 ) .and. ( ftweight(L,ft,iv+1)  <  ftweight(L,ft,iv) ))then
                        !where there is a partly empty leaf layer, some fluxes go straight through.
                        lai_change(L,ft,iv) = ftweight(L,ft,iv)-ftweight(L,ft,iv+1)
                     endif
                     if (ftweight(L,ft,iv+1) - ftweight(L,ft,iv) > 1.e-10_r8)then
                        write(fates_log(),*) 'lower layer has more coverage. This is wrong' , &
                             ftweight(L,ft,iv),ftweight(L,ft,iv+1),ftweight(L,ft,iv+1)-ftweight(L,ft,iv)
                     endif
                     
                     !n.b. in theory lai_change could be calculated daily in the ED code.
                     !This is light coming striaght through the canopy.
                     if (L==1)then
                        tr_dir_z(L,ft,iv+1) = exp(-k_dir(ft) * laisum)* &
                             (ftweight(L,ft,iv)/ftweight(L,ft,1))
                     else   
                        tr_dir_z(L,ft,iv+1) = weighted_dir_tr(L-1)*exp(-k_dir(ft) * laisum)* &
                             (ftweight(L,ft,iv)/ftweight(L,ft,1))
                     endif
                     
                     if (iv == 1)then
                        !this is the top layer.
                        tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * &
                             ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1))                               
                     else
                        !the lai_change(iv) affects the light incident on layer iv+2 not iv+1
                        ! light coming from the layer above (iv-1) goes through iv and onto iv+1. 
                        if (lai_change(L,ft,iv-1) > 0.0_r8)then
                           tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv)* &
                                lai_change(L,ft,iv-1) / ftweight(L,ft,1)
                           tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv-1)* &
                                (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1)
                        else
                           !account fot the light that comes striaght down from unfilled layers above.
                           tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * &
                                ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1))
                        endif
                     endif

                  end do
                  
                  !add up all the weighted contributions from the different PFT columns.
                  weighted_dir_tr(L) = weighted_dir_tr(L) + tr_dir_z(L,ft,currentPatch%nrad(L,ft)+1)*ftweight(L,ft,1)

                  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                  ! Sunlit and shaded fraction of leaf layer
                  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                                
                  !laisum = 0._r8
                  do iv = 1,currentPatch%nrad(L,ft)
                     ! Cumulative leaf area. Original code uses cumulative lai do layer.
                     ! Now use cumulative lai at center of layer.
                     ! Same as tr_dir_z calcualtions, but in the middle of the layer? FIX(RF,032414)-WHY?
                     if (iv  ==  1) then
                        laisum = 0.5_r8 * (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv))
                     else
                        laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv)
                     end if
                     
                     
                     if (L == 1)then !top canopy layer
                        currentPatch%f_sun(L,ft,iv) = exp(-k_dir(ft) * laisum)* &
                             (ftweight(L,ft,iv)/ftweight(L,ft,1))
                     else
                        currentPatch%f_sun(L,ft,iv) = weighted_fsun(L-1)* exp(-k_dir(ft) * laisum)* &
                             (ftweight(L,ft,iv)/ftweight(L,ft,1))
                     endif
                     
                     if ( iv > 1 ) then  ! becasue we are looking at this layer (not the next)
                        ! we only ever add fluxes if iv>1
                        if (lai_change(L,ft,iv-1) > 0.0_r8)then
                           currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + &
                                currentPatch%f_sun(L,ft,iv) * &
                                lai_change(L,ft,iv-1)/ftweight(L,ft,1)
                           currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + &
                                currentPatch%f_sun(L,ft,iv-1) * &
                                (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1)
                        else
                           currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + &
                                currentPatch%f_sun(L,ft,iv-1) * &
                                (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)
                        endif
                     endif
                                   
                  end do !iv

                  weighted_fsun(L) = weighted_fsun(L) + currentPatch%f_sun(L,ft,currentPatch%nrad(L,ft))* &
                       ftweight(L,ft,1)
                  
                  ! instance where the first layer ftweight is used a proxy for the whole column. FTWA
                  ! this is possibly a source of slight error. If we use the ftweight at the top of the PFT column,
                  ! then we willl underestimate fsun, but if we use ftweight at the bottom of the column, we will
                  ! underestimate it. Really, we should be tracking the release of direct light from the column as it tapers
                  ! towards the ground. Is that necessary to get energy closure? It would be quite hard...
               endif !present.
            end do!pft loop
         end do !L
         
         
         do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up.
            do ft = 1,numpft
               if (currentPatch%canopy_mask(L,ft) == 1)then

                  !==============================================================================!
                  ! Iterative solution do scattering
                  !==============================================================================!
                  
                  do ib = 1,hlm_numSWb !vis, nir
                     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                     ! Leaf scattering coefficient and terms do diffuse radiation reflected
                     ! and transmitted by a layer
                     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                     f_not_abs(ft,ib) = rhol(ft,ib) + taul(ft,ib) !leaf level fraction NOT absorbed.
                     !tr_dif_z is a term that uses the LAI in each layer, whereas rhol and taul do not,
                     !because they are properties of leaf surfaces and not of the leaf matrix.
                     do iv = 1,currentPatch%nrad(L,ft)
                        !How much diffuse light is intercepted and then reflected?
                        refl_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * rhol(ft,ib)
                        !How much diffuse light in this layer is transmitted?
                        tran_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * &
                             taul(ft,ib) + tr_dif_z(L,ft,iv)
                     end do
                     
                     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                     ! Ratio of upward to forward diffuse fluxes, dif_ratio
                     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                     ! Soil diffuse reflectance (ratio of down to up radiation).
                     iv = currentPatch%nrad(L,ft) + 1
                     if (L  == currentPatch%NCL_p)then !nearest the soil
                        dif_ratio(L,ft,iv,ib) = currentPatch%gnd_alb_dif(ib)  !bc_in(s)%albgr_dif_rb(ib)
                     else
                        dif_ratio(L,ft,iv,ib) = weighted_dif_ratio(L+1,ib)
                     end if
                     ! Canopy layers, working upwardfrom soil with dif_ratio(iv+1) known
                     ! FIX(RF,032414) ray tracing eqution - need to find derivation of this...
                     ! for each unit going down, there are x units going up.
                     do iv = currentPatch%nrad(L,ft),1, -1
                        dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv+1,ib) * &
                             tran_dif(L,ft,iv,ib)*tran_dif(L,ft,iv,ib) / &
                             (1._r8 - dif_ratio(L,ft,iv+1,ib) * refl_dif(L,ft,iv,ib)) &
                             + refl_dif(L,ft,iv,ib)
                        dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) * &
                             ftweight(L,ft,iv)/ftweight(L,ft,1)
                        dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) + dif_ratio(L,ft,iv+1,ib) * &
                             (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)
                     end do
                     weighted_dif_ratio(L,ib) = weighted_dif_ratio(L,ib) + &
                          dif_ratio(L,ft,1,ib) * ftweight(L,ft,1)
                     !instance where the first layer ftweight is used a proxy for the whole column. FTWA
                  end do!hlm_numSWb
               endif ! currentPatch%canopy_mask
            end do!ft
         end do!L
                       
         
         do ib = 1,hlm_numSWb
            Dif_dn(:,:,:) = 0.00_r8
            Dif_up(:,:,:) = 0.00_r8
            do L = 1, currentPatch%NCL_p !work down from the top of the canopy.
               weighted_dif_down(L) = 0._r8
               do ft = 1, numpft
                  if (currentPatch%canopy_mask(L,ft) == 1)then
                     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                     ! First estimates do downward and upward diffuse flux
                     !
                     ! Dif_dn =  forward diffuse flux onto layer J
                     ! Dif_up =  Upward diffuse flux above layer J
                     !
                     ! Solved here without direct beam radiation and using dif_ratio = Dif_up / Dif_dn
                     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
                     ! downward diffuse flux onto the top surface of the canopy
                     
                     if (L == 1)then
                        Dif_dn(L,ft,1) = forc_dif(radtype)
                     else
                        Dif_dn(L,ft,1) = weighted_dif_down(L-1)
                     end if
                     ! forward diffuse flux within the canopy and at soil, working forward through canopy
                     do iv = 1,currentPatch%nrad(L,ft)
                        denom = refl_dif(L,ft,iv,ib) *  dif_ratio(L,ft,iv,ib)
                        denom = 1._r8 - denom
                        Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) / &
                             denom *ftweight(L,ft,iv)/ftweight(L,ft,1)
                        if (iv > 1)then
                           if (lai_change(L,ft,iv-1) > 0.0_r8)then
                              !here we are thinking about whether the layer above had an laichange,
                              !but calculating the flux onto the layer below.
                              Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv)* &
                                   lai_change(L,ft,iv-1)/ftweight(L,ft,1)
                              Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv-1)* &
                                   (ftweight(L,ft,1)-ftweight(L,ft,iv-1)/ftweight(L,ft,1))
                           else
                              Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * &
                                   (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)
                           endif
                        else
                           Dif_dn(L,ft,iv+1)    = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * &
                                (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)
                        endif
                     end do
                     
                     weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * &
                          ftweight(L,ft,1)
                                   
                     !instance where the first layer ftweight is used a proxy for the whole column. FTWA
                  endif !present
               end do !ft
               if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is the the (incomplete) understorey?
                  !Add on the radiation going through the canopy gaps.
                  weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1)))
                  !instance where the first layer ftweight is used a proxy for the whole column. FTWA
               endif
            end do !L
            
            do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom.
               weighted_dif_up(L) = 0._r8
               do ft = 1, numpft
                  if (currentPatch%canopy_mask(L,ft) == 1)then
                     !Bounce diffuse radiation off soil surface.
                     iv = currentPatch%nrad(L,ft) + 1
                     if (L==currentPatch%NCL_p)then !is this the bottom layer ?
                        Dif_up(L,ft,iv) = currentPatch%gnd_alb_dif(ib) * Dif_dn(L,ft,iv)
                     else
                        Dif_up(L,ft,iv) = weighted_dif_up(L+1)
                     end if
                     ! Upward diffuse flux within the canopy and above the canopy, working upward through canopy
                     
                     do iv = currentPatch%nrad(L,ft), 1, -1
                        if (lai_change(L,ft,iv) > 0.0_r8)then
                           Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * &
                                ftweight(L,ft,iv) / ftweight(L,ft,1)
                           Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * &
                                tran_dif(L,ft,iv,ib) * lai_change(L,ft,iv)/ftweight(L,ft,1)
                           Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * &
                                (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)
                           !nb is this the right constuction?
                           ! the radiation that hits the empty space is not reflected.
                        else
                           Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * ftweight(L,ft,iv)
                           Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * (1.0_r8-ftweight(L,ft,iv))
                        endif
                     end do
                     
                     weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1)
                     !instance where the first layer ftweight is used a proxy for the whole column. FTWA
                  endif !present
               end do !ft
               if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey?
                  !Add on the radiation coming up through the canopy gaps.
                  !diffuse to diffuse
                  weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * &
                       weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib) 
                  !direct to diffuse
                  weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(radtype) * &
                       weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) * currentPatch%gnd_alb_dir(ib)
               endif
            end do !L
            
            !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
            ! 3. Iterative calculation of forward and upward diffuse fluxes, iNCL_puding
            ! scattered direct beam
            !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

            ! Flag to exit iteration loop: 0 = exit and 1 = iterate
            irep = 1
            ! Iteration loop
            iter = 0
            do while(irep ==1 .and. iter<50)
               
               iter = iter + 1
               irep = 0
               do L = 1,currentPatch%NCL_p !working from the top down
                  weighted_dif_down(L) = 0._r8
                  do ft =1,numpft
                     if (currentPatch%canopy_mask(L,ft) == 1)then
                        ! forward diffuse flux within the canopy and at soil, working forward through canopy
                        ! with Dif_up -from previous iteration-. Dif_dn(1) is the forward diffuse flux onto the canopy.
                        ! Note: down = forward flux onto next layer
                        if (L == 1)then !is this the top layer?
                           Dif_dn(L,ft,1) = forc_dif(radtype)
                        else
                           Dif_dn(L,ft,1) = weighted_dif_down(L-1)
                        end if
                        down_rad = 0._r8
                        
                        do iv = 1, currentPatch%nrad(L,ft)
                           
                           down_rad = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) + &
                                Dif_up(L,ft,iv+1) * refl_dif(L,ft,iv,ib) + &
                                forc_dir(radtype) * tr_dir_z(L,ft,iv) * (1.00_r8 - &
                                exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ &
                                currentPatch%esai_profile(L,ft,iv)))) * taul(ft,ib)
                           down_rad = down_rad *(ftweight(L,ft,iv)/ftweight(L,ft,1))
                                         
                           if (iv > 1)then
                              if (lai_change(L,ft,iv-1) > 0.0_r8)then
                                 down_rad = down_rad + Dif_dn(L,ft,iv)   * lai_change(L,ft,iv-1)/ftweight(L,ft,1)
                                 down_rad = down_rad + Dif_dn(L,ft,iv-1) * (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ &
                                      ftweight(L,ft,1)
                              else
                                 down_rad = down_rad + Dif_dn(L,ft,iv)   * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ &
                                      ftweight(L,ft,1)
                              endif
                           else
                              down_rad = down_rad + Dif_dn(L,ft,iv)   * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ &
                                   ftweight(L,ft,1)
                           endif
                                         
                           !this is just Dif down, plus refl up, plus dir intercepted and turned into dif... ,
                           if (abs(down_rad - Dif_dn(L,ft,iv+1)) > tolerance)then
                              irep = 1
                           end if
                           Dif_dn(L,ft,iv+1) = down_rad
                           
                        end do !iv
                        
                        weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * &
                             ftweight(L,ft,1)
                        
                     endif !present
                  end do!ft
                  if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey?
                     weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1) * &
                          (1.0_r8-sum(ftweight(L,1:numpft,1)))
                  end if
               end do ! do L loop

               do L = 1, currentPatch%NCL_p ! working from the top down.
                  weighted_dif_up(L) = 0._r8
                  do ft =1,numpft
                     if (currentPatch%canopy_mask(L,ft) == 1)then
                        ! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam)
                        iv = currentPatch%nrad(L,ft) + 1
                        if (L==currentPatch%NCL_p)then  !In the bottom canopy layer, reflect off the soil
                           Dif_up(L,ft,iv) = Dif_dn(L,ft,iv) * currentPatch%gnd_alb_dif(ib) + &
                                forc_dir(radtype) * tr_dir_z(L,ft,iv) * currentPatch%gnd_alb_dir(ib) 
                        else      !In the other canopy layers, reflect off the underlying vegetation.
                           Dif_up(L,ft,iv) =  weighted_dif_up(L+1)
                        end if
                        
                        ! Upward diffuse flux within and above the canopy, working upward through canopy
                        ! with Dif_dn from previous interation.  Note: up = upward flux above current layer
                        do iv = currentPatch%nrad(L,ft),1,-1
                           !this is radiation up, by layer transmittance, by
                           
                           !reflection of the lower layer,
                           up_rad = Dif_dn(L,ft,iv) * refl_dif(L,ft,iv,ib)
                           up_rad = up_rad + forc_dir(radtype) * tr_dir_z(L,ft,iv) * (1.00_r8 - exp(-k_dir(ft) * &
                                (currentPatch%elai_profile(L,ft,iv) + currentPatch%esai_profile(L,ft,iv)))) * &
                                rhol(ft,ib)
                           up_rad = up_rad + Dif_up(L,ft,iv+1) * tran_dif(L,ft,iv,ib)
                           up_rad = up_rad * ftweight(L,ft,iv)/ftweight(L,ft,1)
                           up_rad = up_rad + Dif_up(L,ft,iv+1) *(ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)
                           ! THE LOWER LAYER FLUX IS HOMOGENIZED, SO WE DON"T CONSIDER THE LAI_CHANGE HERE...
                           
                           if (abs(up_rad - Dif_up(L,ft,iv)) > tolerance) then !are we close to the tolerance level?
                              irep = 1
                           end if
                           Dif_up(L,ft,iv) = up_rad
                           
                        end do  !iv
                        weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1)
                     end if !present
                  end do!ft
                  
                  if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then  !is this the (incomplete) understorey?
                     !Add on the radiation coming up through the canopy gaps.
                     weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * &
                          weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib)
                     weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(radtype) * &
                          weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))*currentPatch%gnd_alb_dir(ib)
                  end if
               end do!L
            end do ! do while over iter
                          
            abs_rad(ib) = 0._r8
            tr_soili = 0._r8
            tr_soild = 0._r8
            
            do L = 1, currentPatch%NCL_p !working from the top down.
               abs_dir_z(:,:) = 0._r8
               abs_dif_z(:,:) = 0._r8
               do ft =1,numpft
                  if (currentPatch%canopy_mask(L,ft) == 1)then
                     !==============================================================================!
                     ! Compute absorbed flux densities
                     !==============================================================================!
                     
                     ! Absorbed direct beam and diffuse do leaf layers
                     do iv = 1, currentPatch%nrad(L,ft)
                        Abs_dir_z(ft,iv) = ftweight(L,ft,iv)* forc_dir(radtype) * tr_dir_z(L,ft,iv) * &
                             (1.00_r8 - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ &
                             currentPatch%esai_profile(L,ft,iv)))) * (1.00_r8 - f_not_abs(ft,ib))
                        Abs_dif_z(ft,iv) = ftweight(L,ft,iv)* ((Dif_dn(L,ft,iv) + &
                             Dif_up(L,ft,iv+1)) * (1.00_r8 - tr_dif_z(L,ft,iv)) * &
                             (1.00_r8 - f_not_abs(ft,ib)))
                     end do
                     
                     ! Absorbed direct beam and diffuse do soil
                     if (L == currentPatch%NCL_p)then
                        iv = currentPatch%nrad(L,ft) + 1
                        Abs_dif_z(ft,iv) = ftweight(L,ft,1)*Dif_dn(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dif(ib) )
                        Abs_dir_z(ft,iv) = ftweight(L,ft,1)*forc_dir(radtype) * &
                             tr_dir_z(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dir(ib)  )
                        tr_soild = tr_soild + ftweight(L,ft,1)*forc_dir(radtype) * tr_dir_z(L,ft,iv)
                        tr_soili = tr_soili + ftweight(L,ft,1)*Dif_dn(L,ft,iv)
                     end if

                     ! Absorbed radiation, shaded and sunlit portions of leaf layers
                     !here we get one unit of diffuse radiation... how much of
                     !it is absorbed?
                     if (ib == ivis) then ! only set the absorbed PAR for the visible light band. 
                        do iv = 1, currentPatch%nrad(L,ft)
                           if (radtype==idirect) then
                              if ( debug ) then
                                 write(fates_log(),*) 'EDsurfAlb 730 ',Abs_dif_z(ft,iv),currentPatch%f_sun(L,ft,iv)
                                 write(fates_log(),*) 'EDsurfAlb 731 ', currentPatch%fabd_sha_z(L,ft,iv), &
                                      currentPatch%fabd_sun_z(L,ft,iv)
                              endif
                              currentPatch%fabd_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * &
                                   (1._r8 - currentPatch%f_sun(L,ft,iv))
                              currentPatch%fabd_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * &
                                   currentPatch%f_sun(L,ft,iv) + &
                                   Abs_dir_z(ft,iv)
                           else
                              currentPatch%fabi_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * &
                                   (1._r8 - currentPatch%f_sun(L,ft,iv))
                              currentPatch%fabi_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * &
                                   currentPatch%f_sun(L,ft,iv)
                           endif
                           if ( debug ) then
                              write(fates_log(),*) 'EDsurfAlb 740 ', currentPatch%fabd_sha_z(L,ft,iv), &
                                   currentPatch%fabd_sun_z(L,ft,iv)
                           endif
                        end do
                     endif ! ib 

                     
                     !==============================================================================!
                     ! Sum fluxes
                     !==============================================================================!
                     ! Solar radiation absorbed by ground
                     iv = currentPatch%nrad(L,ft) + 1
                     if (L==currentPatch%NCL_p)then
                        abs_rad(ib) = abs_rad(ib) +  (Abs_dir_z(ft,iv) + Abs_dif_z(ft,iv))
                     end if
                     ! Solar radiation absorbed by vegetation and sunlit/shaded leaves
                     do iv = 1,currentPatch%nrad(L,ft)
                        if (radtype == idirect)then
                           currentPatch%fabd(ib) = currentPatch%fabd(ib) + &
                                Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv)
                           ! bc_out(s)%fabd_parb_out(ib) = currentPatch%fabd(ib)
                        else
                           currentPatch%fabi(ib) = currentPatch%fabi(ib) + Abs_dif_z(ft,iv)
                           ! bc_out(s)%fabi_parb_out(ib) = currentPatch%fabi(ib)
                        endif
                     end do
                     
                     ! Albefor
                     if (L==1)then !top canopy layer.
                        if (radtype == idirect)then
                           albd_parb_out(ib) = albd_parb_out(ib) + &
                                Dif_up(L,ft,1) * ftweight(L,ft,1)
                        else
                           albi_parb_out(ib) = albi_parb_out(ib) + &
                                Dif_up(L,ft,1) * ftweight(L,ft,1)
                        end if
                     end if
                     
                     ! pass normalized PAR profiles for use in diagnostic averaging for history fields
                     if (ib == ivis) then ! only diagnose PAR profiles for the visible band
                        do iv = 1, currentPatch%nrad(L,ft)
                           currentPatch%nrmlzd_parprof_pft_dir_z(radtype,L,ft,iv) = &
                                forc_dir(radtype) * tr_dir_z(L,ft,iv)
                           currentPatch%nrmlzd_parprof_pft_dif_z(radtype,L,ft,iv) = &
                                Dif_dn(L,ft,iv) + Dif_up(L,ft,iv)
                           !
                           currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) = &
                                currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) + &
                                (forc_dir(radtype) * tr_dir_z(L,ft,iv)) * &
                                (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv)))
                           currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) = &
                                currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) + &
                                (Dif_dn(L,ft,iv) + Dif_up(L,ft,iv)) * &
                                (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv)))
                        end do
                     end if ! ib = visible
                  end if ! present
               end do !ft
               if (radtype == idirect)then
                  fabd_parb_out(ib) = currentPatch%fabd(ib)
               else
                  fabi_parb_out(ib) = currentPatch%fabi(ib)
               endif
               
               
               !radiation absorbed from fluxes through unfilled part of lower canopy.
               if (currentPatch%NCL_p > 1.and.L == currentPatch%NCL_p)then 
                  abs_rad(ib) = abs_rad(ib) + weighted_dif_down(L-1) * &
                       (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dif(ib) )
                  abs_rad(ib) = abs_rad(ib) + forc_dir(radtype) * weighted_dir_tr(L-1) * &
                       (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dir(ib) )
                  tr_soili = tr_soili + weighted_dif_down(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))
                  tr_soild = tr_soild + forc_dir(radtype) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))
               endif
               
               if (radtype == idirect)then
                  currentPatch%tr_soil_dir(ib) = tr_soild
                  currentPatch%tr_soil_dir_dif(ib) = tr_soili
                  currentPatch%sabs_dir(ib)     = abs_rad(ib)
                  ftdd_parb_out(ib)  = tr_soild
                  ftid_parb_out(ib) =  tr_soili
               else
                  currentPatch%tr_soil_dif(ib) = tr_soili
                  currentPatch%sabs_dif(ib)     = abs_rad(ib)
                  ftii_parb_out(ib) =  tr_soili
               end if
               
            end do!l
            
                          
            !==============================================================================!
            ! Conservation check
            !==============================================================================!
            ! Total radiation balance: absorbed = incoming - outgoing
            
            if (radtype == idirect)then
               error = abs(currentPatch%sabs_dir(ib) - (currentPatch%tr_soil_dir(ib) * &
                    (1.0_r8-currentPatch%gnd_alb_dir(ib) ) + &
                    currentPatch%tr_soil_dir_dif(ib) * (1.0_r8-currentPatch%gnd_alb_dif(ib)     )))
               if ( abs(error) > 0.0001)then
                  write(fates_log(),*)'dir ground absorption error',error,currentPatch%sabs_dir(ib), &
                       currentPatch%tr_soil_dir(ib)* &
                       (1.0_r8-currentPatch%gnd_alb_dir(ib)  ),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1))
                  write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), &
                       (1.0_r8-currentPatch%gnd_alb_dir(ib)  )
                  
                  do ft =1,3
                     iv = currentPatch%nrad(1,ft) + 1
                     write(fates_log(),*) 'abs soil fluxes', Abs_dir_z(ft,iv),Abs_dif_z(ft,iv)
                  end do
                  
               end if
            else
               if ( abs(currentPatch%sabs_dif(ib)-(currentPatch%tr_soil_dif(ib) * &
                    (1.0_r8-currentPatch%gnd_alb_dif(ib)  ))) > 0.0001_r8)then
                  write(fates_log(),*)'dif ground absorption error',currentPatch%sabs_dif(ib) , &
                       (currentPatch%tr_soil_dif(ib)* &
                       (1.0_r8-currentPatch%gnd_alb_dif(ib)  )),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1))
               endif
            endif
            
            if (radtype == idirect)then
               error = (forc_dir(radtype) + forc_dif(radtype)) - &
                    (fabd_parb_out(ib)  + albd_parb_out(ib) + currentPatch%sabs_dir(ib))
            else
               error = (forc_dir(radtype) + forc_dif(radtype)) - &
                    (fabi_parb_out(ib)  + albi_parb_out(ib) + currentPatch%sabs_dif(ib))
            endif
            lai_reduction(:) = 0.0_r8
            do L = 1, currentPatch%NCL_p
               do ft =1,numpft
                  if (currentPatch%canopy_mask(L,ft) == 1)then
                     do iv = 1, currentPatch%nrad(L,ft)
                        if (lai_change(L,ft,iv) > 0.0_r8)then
                           lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv))
                        endif
                     enddo
                  endif
               enddo
            enddo
            
            if (radtype == idirect)then
               !here we are adding a within-ED radiation scheme tolerance, and then adding the diffrence onto the albedo
               !it is important that the lower boundary for this is ~1000 times smaller than the tolerance in surface albedo. 
               if (abs(error)  >  1.e-9_r8 .and. abs(error) < 0.15_r8)then
                  albd_parb_out(ib) = albd_parb_out(ib) + error
                  !this terms adds the error back on to the albedo. While this is partly inexcusable, it is 
                  ! in the medium term a solution that
                  ! prevents the model from crashing with small and occasional energy balances issues.
                  ! These are extremely difficult to debug, many have been solved already, leading
                  ! to the complexity of this code, but where the system generates occasional errors, we
                  ! will deal with them for now.
               end if
               if (abs(error)  >  0.15_r8)then
                  write(fates_log(),*) 'Large Dir Radn consvn error',error ,ib
                  write(fates_log(),*) 'diags', albd_parb_out(ib), ftdd_parb_out(ib), &
                       ftid_parb_out(ib), fabd_parb_out(ib)
                  write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'ftweight',ftweight(1,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno
                  write(fates_log(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib)
                  
                  albd_parb_out(ib) = albd_parb_out(ib) + error
               end if
            else
               
               if (abs(error)  >  1.e-9_r8 .and. abs(error) < 0.15_r8)then
                  albi_parb_out(ib) = albi_parb_out(ib) + error
               end if
               
               if (abs(error)  >  0.15_r8)then
                  write(fates_log(),*)  '>5% Dif Radn consvn error',error ,ib
                  write(fates_log(),*) 'diags', albi_parb_out(ib), ftii_parb_out(ib), &
                       fabi_parb_out(ib)
                  write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf)
                  write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno
                  write(fates_log(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib)
                  write(fates_log(),*) 'rhol',rhol(1:numpft,:)
                  write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1)
                  write(fates_log(),*) 'present',currentPatch%canopy_mask(1,1:numpft)
                  write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1)
                  
                  albi_parb_out(ib) = albi_parb_out(ib) + error
               end if
               
               if (radtype == idirect)then
                  error = (forc_dir(radtype) + forc_dif(radtype)) - &
                       (fabd_parb_out(ib)  + albd_parb_out(ib) + currentPatch%sabs_dir(ib))
               else
                  error = (forc_dir(radtype) + forc_dif(radtype)) - &
                       (fabi_parb_out(ib)  + albi_parb_out(ib) + currentPatch%sabs_dif(ib))
               endif
               
               if (abs(error)  >  0.00000001_r8)then
                  write(fates_log(),*)  'there is still error after correction',error ,ib
               end if
               
            end if
            
         end do !hlm_numSWb
         
      enddo ! rad-type
      
      
    end associate
    return
  end subroutine PatchNormanRadiation

 ! ======================================================================================

 subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out)
    
    implicit none

    ! Arguments
    integer,intent(in)                      :: nsites
    type(ed_site_type),intent(inout),target :: sites(nsites)
    type(bc_in_type),intent(in)             :: bc_in(nsites)
    type(bc_out_type),intent(inout)         :: bc_out(nsites)
    

    ! locals
    type (ed_patch_type),pointer :: cpatch   ! c"urrent" patch
    real(r8)          :: sunlai
    real(r8)          :: shalai
    real(r8)          :: elai
    integer           :: CL
    integer           :: FT
    integer           :: iv
    integer           :: s
    integer           :: ifp
    

    do s = 1,nsites

       ifp = 0
       cpatch => sites(s)%oldest_patch

       do while (associated(cpatch))                 
          
          ifp=ifp+1
          
          if( debug ) write(fates_log(),*) 'edsurfRad_5600',ifp,s,cpatch%NCL_p,numpft
          
          ! zero out various datas
          cpatch%ed_parsun_z(:,:,:) = 0._r8
          cpatch%ed_parsha_z(:,:,:) = 0._r8
          cpatch%ed_laisun_z(:,:,:) = 0._r8     
          cpatch%ed_laisha_z(:,:,:) = 0._r8

          bc_out(s)%fsun_pa(ifp) = 0._r8

          sunlai  = 0._r8
          shalai  = 0._r8

          cpatch%parprof_pft_dir_z(:,:,:) = 0._r8
          cpatch%parprof_pft_dif_z(:,:,:) = 0._r8
          cpatch%parprof_dir_z(:,:) = 0._r8
          cpatch%parprof_dif_z(:,:) = 0._r8

          ! Loop over patches to calculate laisun_z and laisha_z for each layer.
          ! Derive canopy laisun, laisha, and fsun from layer sums.
          ! If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from
          ! SurfaceAlbedo is canopy integrated so that layer value equals canopy value.
          
          ! cpatch%f_sun is calculated in the surface_albedo routine...
          
          do CL = 1, cpatch%NCL_p
             do FT = 1,numpft

                if( debug ) write(fates_log(),*) 'edsurfRad_5601',CL,FT,cpatch%nrad(CL,ft)
                
                do iv = 1, cpatch%nrad(CL,ft) !NORMAL CASE. 
                   
                   ! FIX(SPM,040114) - existing comment
                   ! ** Should this be elai or tlai? Surely we only do radiation for elai? 
                   
                   cpatch%ed_laisun_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * &
                         cpatch%f_sun(CL,ft,iv)
                   
                   if ( debug ) write(fates_log(),*) 'edsurfRad 570 ',cpatch%elai_profile(CL,ft,iv)
                   if ( debug ) write(fates_log(),*) 'edsurfRad 571 ',cpatch%f_sun(CL,ft,iv)
                   
                   cpatch%ed_laisha_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * &
                         (1._r8 - cpatch%f_sun(CL,ft,iv))
                   
                end do
                
                !needed for the VOC emissions, etc. 
                sunlai = sunlai + sum(cpatch%ed_laisun_z(CL,ft,1:cpatch%nrad(CL,ft)))
                shalai = shalai + sum(cpatch%ed_laisha_z(CL,ft,1:cpatch%nrad(CL,ft)))
                
             end do
          end do
          
          if(sunlai+shalai > 0._r8)then
             bc_out(s)%fsun_pa(ifp) = sunlai / (sunlai+shalai) 
          else
             bc_out(s)%fsun_pa(ifp) = 0._r8
          endif
          
          if(bc_out(s)%fsun_pa(ifp) > 1._r8)then
             write(fates_log(),*) 'too much leaf area in profile',  bc_out(s)%fsun_pa(ifp), &
                   sunlai,shalai
          endif

          elai = calc_areaindex(cpatch,'elai')

          bc_out(s)%laisun_pa(ifp) = elai*bc_out(s)%fsun_pa(ifp)
          bc_out(s)%laisha_pa(ifp) = elai*(1.0_r8-bc_out(s)%fsun_pa(ifp))

         ! Absorbed PAR profile through canopy
         ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo
         ! are canopy integrated so that layer values equal big leaf values.
         
         if ( debug ) write(fates_log(),*) 'edsurfRad 645 ',cpatch%NCL_p,numpft
         
         do CL = 1, cpatch%NCL_p
            do FT = 1,numpft
               
               if ( debug ) write(fates_log(),*) 'edsurfRad 649 ',cpatch%nrad(CL,ft)
               
               do iv = 1, cpatch%nrad(CL,ft)
                  
                  if ( debug ) then
                     write(fates_log(),*) 'edsurfRad 653 ', cpatch%ed_parsun_z(CL,ft,iv)
                     write(fates_log(),*) 'edsurfRad 654 ', bc_in(s)%solad_parb(ifp,ipar)
                     write(fates_log(),*) 'edsurfRad 655 ', bc_in(s)%solai_parb(ifp,ipar)
                     write(fates_log(),*) 'edsurfRad 656 ', cpatch%fabd_sun_z(CL,ft,iv)
                     write(fates_log(),*) 'edsurfRad 657 ', cpatch%fabi_sun_z(CL,ft,iv)
                  endif
                  
                  cpatch%ed_parsun_z(CL,ft,iv) = &
                        bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sun_z(CL,ft,iv) + &
                        bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sun_z(CL,ft,iv) 
                  
                  if ( debug )write(fates_log(),*) 'edsurfRad 663 ', cpatch%ed_parsun_z(CL,ft,iv)
                  
                  cpatch%ed_parsha_z(CL,ft,iv) = &
                        bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sha_z(CL,ft,iv) + &
                        bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sha_z(CL,ft,iv)          
                  
                  if ( debug ) write(fates_log(),*) 'edsurfRad 669 ', cpatch%ed_parsha_z(CL,ft,iv)
                  
               end do !iv
            end do !FT
         end do !CL

         ! output the actual PAR profiles through the canopy for diagnostic purposes
         
         do CL = 1, cpatch%NCL_p
            do FT = 1,numpft
               do iv = 1, cpatch%nrad(CL,ft)
                  cpatch%parprof_pft_dir_z(CL,FT,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
                       cpatch%nrmlzd_parprof_pft_dir_z(idirect,CL,FT,iv)) + &
                       (bc_in(s)%solai_parb(ifp,ipar) * &
                       cpatch%nrmlzd_parprof_pft_dir_z(idiffuse,CL,FT,iv))
                  cpatch%parprof_pft_dif_z(CL,FT,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
                       cpatch%nrmlzd_parprof_pft_dif_z(idirect,CL,FT,iv)) + &
                       (bc_in(s)%solai_parb(ifp,ipar) * &
                       cpatch%nrmlzd_parprof_pft_dif_z(idiffuse,CL,FT,iv))
               end do ! iv
            end do    ! FT
         end do       ! CL

         do CL = 1, cpatch%NCL_p
            do iv = 1, maxval(cpatch%nrad(CL,:))
               cpatch%parprof_dir_z(CL,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
                    cpatch%nrmlzd_parprof_dir_z(idirect,CL,iv)) + &
                    (bc_in(s)%solai_parb(ifp,ipar) * &
                    cpatch%nrmlzd_parprof_dir_z(idiffuse,CL,iv))
               cpatch%parprof_dif_z(CL,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
                    cpatch%nrmlzd_parprof_dif_z(idirect,CL,iv)) + &
                    (bc_in(s)%solai_parb(ifp,ipar) * &
                    cpatch%nrmlzd_parprof_dif_z(idiffuse,CL,iv))
            end do    ! iv
         end do       ! CL
         
         cpatch => cpatch%younger
      enddo
      
      
   enddo
   return
   
end subroutine ED_SunShadeFracs


!      ! MOVE TO THE INTERFACE
!      subroutine ED_CheckSolarBalance(g,filter_nourbanp,num_nourbanp,fsa,fsr,forc_solad,forc_solai)


!         implicit none
!         integer,intent(in),dimension(:)    :: gridcell     ! =>    gridcell index
!         integer,intent(in),dimension(:)    :: filter_nourbanp ! => patch filter for non-urban points
!         integer, intent(in)                :: num_nourbanp !       number of patches in non-urban points in patch  filter
!         real(r8),intent(in),dimension(:,:) :: forc_solad   ! =>    atm2lnd_inst%forc_solad_grc, direct radiation (W/m**2
!         real(r8),intent(in),dimension(:,:) :: forc_solai   ! =>    atm2lnd_inst%forc_solai_grc, diffuse radiation (W/m**2)
!         real(r8),intent(in),dimension(:,:) :: fsa          ! =>    solarabs_inst%fsa_patch, solar radiation absorbed (total) (W/m**2)
!         real(r8),intent(in),dimension(:,:) :: fsr          ! =>    solarabs_inst%fsr_patch, solar radiation reflected (W/m**2)      

!         integer :: p
!         integer :: fp
!         integer :: g
!         real(r8) :: errsol

!         do fp = 1,num_nourbanp
!            p = filter_nourbanp(fp)
!            g = gridcell(p)
!            errsol = (fsa(p) + fsr(p)  - (forc_solad(g,1) + forc_solad(g,2) + forc_solai(g,1) + forc_solai(g,2)))
!            if(abs(errsol) > 0.1_r8)then
!               write(fates_log(),*) 'sol error in surf rad',p,g, errsol
!            endif
!         end do
!         return
!      end subroutine ED_CheckSolarBalance
   

end module EDSurfaceRadiationMod
